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1. Introduction and motivation 

A central problem in quantitative finance is the development of efficient methods for 
pricing and hedging derivative securities [U 12 E] ■ Although the classical Black, Scholes 
and Merton model of financial derivatives 0| IH] provides an elegant framework to price 
financial derivatives, its actual analytical tractability is limited to plain vanilla call 
and put options and to few other cases. Actually, even if some particular payoffs lead 
to exact or approximated closed-form pricing formulae [HI Ej , these analytical results 
can be extended to more general payoffs only with difficulty. Hence, there is a need 
for flexible and fast pricing algorithms, especially when we are interested in pricing 
options whose payoff at the expiry date depends on the whole path followed by the 
underlying (i.e. path dependent options). In the past years many approaches have been 
proposed and the standard numerical procedures adopted in financial engineering involve 
the use of binomial or trinomial trees, Monte Carlo simulations and finite difference 
methods PQEHSj- Alternative and more recent algorithms are described, for example, 
in [SI, which has a comprehensive bibliography. 

In this paper we extend the path integral approach to option pricing developed for 
unidimensional assets in [Hj. We generalize the original formulation in order to price a 
variety of commonly traded exotic options. First, we obtain a pricing formula for path 
dependent options based on multiple correlated underlying assets; second, we improve 
the related numerical algorithms. Comparisons with standard Monte Carlo simulations, 
as well as with the results of other numerical techniques known in the literature, are 
presented. Related attempts to price options and, in particular, exotic options, using 
the path integral method can be found in[TnilIIllIllIllIllI3linilIIl- 

The structure of the article is the following. In Section |21 we trace the derivation 
of the central pricing formula of our path integral-inspired approach and we describe in 
Section ini how to implement it numerically. Details about the computational algorithms 
used for pricing and numerical results are discussed in Section |3] In this latter, we show 
that our approach can be efficiently implemented to price a large class of exotic options: 
Asian, barrier knock out and reverse cliquet. Finally, in Section El we compute the 
Greeks relative to some of the considered options and in Section IHl we draw conclusions 
and consider possible perspectives. 

2. A path integral-based pricing formula 

Path integral techniques, widely used in quantum mechanics and quantum field theory, 
can be useful to describe the dynamics of a Markov stochastic process [IHl dl 120] • In 
the present paper we are interested in computing mean values of functionals of a D- 
dimensional Markov stochastic process S{t) corresponding to the price of a set of D 
underlying assets. 

In particular, let us fix a time horizon T > and split the time interval [0, T] into 
n + 1 subintervals [Tj, Tj+i] with Tq = 0, T„+i = T and Tj+i — Tj = At = T/{n+l). Our 
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aim is to compute the fair price at time Tq of a European path dependent option with 
maturity T and whose payoff / is a function of the values S{0), S(Ti), . . . , S{Tn+i). 

According to the arbitrage-free pricing theory jH El E] , this means that we have 
to evaluate the mathematical expectation^ 

E[f{Z{Tn+i),Z{Tn),...,Zo)] = 

(1) 

jRDx(n + l) 

where Z(t) = log5'(t) and p{zn+i, Zn, . . . , zq) is the joint probability density function 
(pdf) of the path {Z{To) = Zo,...,Z{Ti) = Zi, . . . , Z{Tn+i) = Zn+i}. In order to 
compute Equation (0), a straightforward application of standard Monte Carlo estimation 
theory, would require the sampling of independent and identically distributed (i.i.d.) 
paths {Zo, Z(')(Ti), . . . , Z'^^\Tn+i)}i=i,...,N ■ Whenever S{t) (and hence Z{t)) is solution 
of a stochastic differential equation (SDE), this can be done by an Euler discretization 
scheme of the SDE. Unfortunately, this procedure may be slow, as to compute Z^^\Ti) 
we typically need to know Z*^')(Tj_i), and is not efficient when considering out-of-the- 
money (OTM) options, because a relevant number of sampled paths may not contribute 
to the payoff. That is why we look for an alternative formulation of this pricing problem 
leading to a reduction of the Monte Carlo variance"*". 



2.1. The model 

Let us first of all introduce the evolution model we will focus on throughout the rest of 
the paper. We assume that S{t) satisfies, under the objective probability measure, the 
following SDE 

dS^{t)/S^{t) = pi^dt + akdW^it) VA; = !,...,£) 

{dW\t), dW^it)) = pijdt Vz, J = 1, . . . , D, ^ ^ 

where the /x'^ G M and the ak G represent the mean returns and the volatilities of S^, 
respectively, and the pij G (—1,1) are the correlations between the components of the 
Wiener processes W{t) {pa = 1). This formulation is particularly useful because it only 
involves financial quantities that can be historically estimated. For instance, pij and 
a'^ can be evaluated by analyzing the time series of the correlations between different 
assets' returns, i.e. 

{dS\t),dS^{t)) = {S\t)a,fdt 
{dS\t),dS^{t)) = S\t)S^{t)a,(jjPijdt i^j. 

However, it is convenient to write Equation in terms of the square root S of 
the variance-covariance matrix Ej = a^ajPij and of a standard D-dimensional Wiener 
process W . The square root S is defined by relation SS"^ = S and can be chosen 
to be a lower triangular matrix. Changing the probability measure from the original 

^ For simplicity, we have included the discount factor exp{— rT} in the definition of the payofi'. 
"•" We refer to 01221 fo'" ^ review of standard Monte Carlo variance reduction techniques. 
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objective measure of Equation to the risk neutral one the stochastic process 
Z(t) = (log5'i(t), . . . , log Soit)) satisfies the following equation 

dZ{t) = Adt + i:dW{t) 

Z{0) = Zo, ^ ^ 

where the fc*^ entry of A is A'' = (r — 0-^/2), with r the risk-free interest rate. From 
Equation we infer that Z{t) is normally distributed with mean Zo + At and variance- 
covariance matrix St. Equivalently, the conditional pdf * p{z' ,t'\z,t), t' > t, is given 
by 

where || ■ || stands for the standard Euclidean norm. Solutions of Equation are 
Markov processes and, therefore, it is possible to describe their time evolution via a 
path integral formulation 



2.2. The fundamental pricing formula 



Thanks to the properties of the chosen model, we are now able to extend the pricing 
formula given in ^ and propose improved algorithms to evaluate Equation (P). Our 
approach is essentially based on a sequence of linear changes of the integration variables 
appearing in Equation (^. This latter expression will then be rewritten in terms of a 
suitable set of independent random variables. 

The definition of conditional probability, together with the Markov nature of the 
price dynamics, allows us to write the joint probability entering Equation (jT]) as 



n+1 



n+1 



P(2^n+1; ^n; • • • ; 

Zq) W dx 



Y\<^Zi Pi_i{Zi\Zi_i) 
i=l 
n+1 



27rAt 



D/2 



IdetSl 



,-||S-i[^,-(^,_i+AAt)]||2/2At 



(6) 



where we have written the transition densities pi^i{zi\zi^i) = p{zi,Ti\zi-i,Ti^i) 
explicitly. In order to get rid of the correlation matrix S and of the drift A, we perform 
a first change of variable by setting Zi = S(?7j + AiAt), i = to n + 1, thus obtaining 
for the r.h.s. of Equation (jH} 



n+1 



:.h.s. © = n 



1 



D/2 



drji exp 



1 



\Vi -Vi-i\ 



(7) 



2nAt 2At' 

i=i L ' 

* By definition, p(z\t'\z,t) is such that the probability for Z(t') taking a value in the _D-dimensional 
hyper-cube dz' centred on z', conditional on Z{t) — z, is p{z' ,t'\z,t)dz' . 
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We work out the quadratic form X]r=i ll'^* ~ rearrange terms by 

introducing the D-dimensional vectors hi, . . . ,hn such that 



= k = l,...D, i = l,...,n, 

where O is the orthogonal matrix that diagonahzes the n x n tridiagonal matrix 

/2 -10 0\ 

-1 2 -10 ■■■0 

-1 2 -1 ■•■ 

... -1 2 -10 

-12-1 

VO -12/ 
After some tedious, but straightforward, algebra, we obtain 



M 



(8) 



r.h.s. (0) = g{zn+i; zo)dzn+i JJ <lhiQi{hu ^o, 



Zn+l] 



(9) 



i=l 



where the Qi{ ■ ; Zq, 2„+i) are D-dimensional Gaussian pdfs with mean di = [H^^zqOu + 

T,~^{zn+i — AiAt)Oni]/mi and variance {At/mi)l[r)xD, the {mj}j=i_...^„ being the 

eigenvalues of M. The function g is defined as 

/ \ D 

^ I 1 \ 

g{Zn+l'-, Zq) 



|detS| \^27rAMet(M) 

|S-i^o|P + - A(n + l)At)||2 - Er=i ll^dl'M 



exp 



2At 



and from now on we will drop, for simplicity, its dependence on zq. 

Finally, we replace the hi by setting hi = 'di + Xi^/ At/mi, thus obtaining the 
ultimate relationship 



r.h.s. (0) = g{zn+i)dzn+iY\_(^^iPG{K 



(10) 



1=1 



where pc is a D-dimensional standardized Gaussian pdf. By means of this sequence of 
replacements, the expectation ((T)) can be computed as 



E[/(Z(T„+i),...,Zo)] = 

(^Zn-\-ig{Zn+l] 



n. n 

I TT (dAiPG(Ai)) f{zn+i, A„, . . . , Ai, zo) 



where f{zn+i, A„, . . . , Ai, zq) = f{zn+i, Zn, . . . , zi, zq) with the substitution 

/At 



iAAt, i = 1 n. 



(11) 



(12) 



z, = Y,Oi,T.[\j—\+^, 

i=i 

The reformulation of Equation (0) as in Equation (fTT|) is the core of our pricing 
technique. In particular, advantages come from having split the D x (n + l)-dimensional 
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integral into an external integration over Zn+i, representing the value of the log-price 
at the maturity, and an internal one, which can be thought as the mathematical 
expectation 

» n 

E[f{zn+i,An, • • • , Ai, ^o)] = / TT (dXiPciXi)) f{zn+i, A„, . . . , Ai, zq), 

where Ai, . . . , A„ are standardized D-dimensional i.i.d. Gaussian variables. 

Let us stress that, for each value of Zn+i, and since zq is known, by means of 
i.i.d. Gaussian samples < (A,''', . . . , An"*) \ , we can construct the set of log-price 

I J 1=1,.. .,N 

paths 

Z«(T.) = J20i,^ {'^^^^^^^^ ^^^^ 

having fixed starting and end points Z('^(T„_|_i) = Zn+i and Z(')(To) = Zq. This way of 
proceeding is typical of path integral techniques, in which functional trajectories with 
fixed initial and final states are considered [IHl CHI EOj • That is why we call our method 
path integral pricing. 



3. Computational algorithms 

The reformulation of the pricing problem given in Section |21 is useful to price path 
dependent options: this task has been reduced to the numerical computation of the 
integrals appearing in Equation (fTTj) . In particular, we can adopt the two following 
procedures 

1 . We can compute the internal D x n-dimensional integral via Monte Carlo sampling 
of the Aj, and the external /^-dimensional one by a deterministic method to be 
specified. We will call this method path integral with external integration. This 
method turns out to perform well when D = 1, as shown in the following. 

2. We can perform a pure D x (?T, + l)-dimensional Monte Carlo integration by properly 
truncating the integration domain on Zn+i- This method will be called pure Monte 
Carlo and is particularly useful when considering OTM options on multidimensional 
assets. 

We provide, in the next two subsections, a more exhaustive insight into the 
procedures sketched above. We refer to Section 0] for the implemented versions' details 
and the numerical results. 



3.1. Path integral with external integration 

This method corresponds to a very precise evaluation of the inner function 
E[/(z„+i, A„, . . . , Ai, 2;o)] = S{zn+i) for some given values of Zn+i- Actually, our aim 
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is to approximate Equation (jllj) by a formula like 



/ dZn+l g{Zn+l)S{Zn+l) ^ V (14) 

with a suitable choice of the integration weights Wi and of the integration points 2:^+1. 



We can, for example, perform Riemann integration or exploit a quadrature rule 
Since £^(z^*|^) is a non-explicitly solvable mathematical expectation, for each we 
estimate it by sampling N i.i.d. from the law of (Ai, . . . , A„), thus obtaining, at the 
same time, the associated errors fj. By virtue of the Central Limit Theorem, the Vi 
scale with the square root of A^, so that the bigger A^ is, the smaller the error and more 
precise are the values of 8{z^^j^]). Of course, the choice of the z)^^^ influences the final 
result and has to be done carefully. By means of these coupled outer-deterministic and 
inner-Monte Carlo integrations, we estimate the price with 



i3± = ^u.,9(4'l,)f(4'li)± 



(15) 



as boundary values for the 68% confidence interval of Equation |TT| . It is worth noticing 
that such an error does not include the effect of finiteness of Uint- Numerical results 
providing us with evidence the error due to a finite Ui^t is negligible are reported in 
Section 

To conclude, let us remark that this procedure, forcing Z{Tn+i) to take a value 
in . . . ,-2^+^}, is similar to the variance reduction technique known as stratified 

sampling Monte Carlo 



3.2. Pure Monte Carlo 

We will show in the next section that when pricing unidimensional assets according 
to Equation ()15|). a deterministic choice of final integration points works better than 
a Monte Carlo one. However, it is known that deterministic integration approaches 
rapidly lose their competitiveness as the dimension grows. As an alternative, we propose 
a method based on a pure Monte Carlo integration coupled with the path integral. 
First, we choose a function V : (0, +00) such that F > and T{z)dz = 1 and 

we interpret it as a pdf. Second, we rewrite Equation (fTTj) as 

„ n 

E[/(Z(T„+i), . . . , Zq)] = / dzn+iT{zn+i) TT pG{Xi)dXif{zn+i, A„, . . . , Ai, 2:0) 

= E[/(Z,A„,...,Ai,zo)], (16) 

where f{zn+i, A„, . . . , Ai, ^o) = g{zn+i)f{zn+i, A„, . . . , Ai, zo)/T{zn+i) and Z is a random 
variable with F as pdf. In other words, we read the pricing formula as the mathematical 
expectation of a function of n + 1 independent variables, namely Z and the Aj. Our 
algorithm evaluates Equation (|T?)|l by a pure Monte Carlo method extracting A^ random 
i.i.d. samples (Z«, A^ , . . . , A'k^)i=i,...,N. 
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This method resembles a standard Monte Carlo simulation of random walks, but 
there are some subtle differences. First of all, in the random walk case one simulates 
each path recursively by sampling n + 1 Gaussian variables without knowing where the 
considered path will end, while here we want to construct paths leading to a given -z^+i- 
Second, we introduce an asymmetry between Zn+i and the Aj in the sense that Zn+i 
plays a crucial role and we give to it the possibility of being sampled by a non-Gaussian 
pdf by changing F. This turns out to be very useful when pricing OTM options and the 
Monte Carlo random walk turns out to be not efficient, as shown in the next section. 

4. Numerical results and discussion 

In this section we report computational issues concerning the pricing of different kinds 
of path dependent options by means of the path integral procedures discussed in 
Section El We will consider the following types of options: Asian and up-out barrier 
unidimensional call, unidimensional reverse cliquet and Asian basket call. The dynamics 
of the underlying assets is assumed to follow Equation (jH). 

4.I. Algorithms' implemented versions 

Before entering into details with numerical results, let us list here which versions of the 
two general algorithms of Section 01 have been implemented. 

(i) Path Integral with Trapezoidal Integration (PITP) 

This is an algorithm of the type described in Section IH.ll In particular, as 
deterministic method used to integrate over Zn+i, we choose trapezoidal integration 
with equispaced abscissa [23]- Integration, whenever not differently specified, 
instead of being performed on M^, is truncated, for each k = 1 to D, to the 
interval 

Ck = [z" -4akVT,z' + 4akVT], 

where z^ = Zq + {r — a\/2)T for in-the-money (ITM) and at-the-money (ATM) 
options and z = log{K) for OTM ones. Tests regarding this choice can be found in 
Section ion 

(ii) Pure Monte Carlo with Flat Sampling (PIFL) 

This is a Pure Monte Carlo algorithm (see Section in which we sample the Z{1) 
of Equation (jTHl) from a uniform pdf on the compact subset C = Ci ^ ■ ■ ■ ® Cd- 

(iii) Pure Monte Carlo with Truncated Cauchy Sampling (PICH) 

This version of the Pure Monte Carlo method uses, as F, a truncated Cauchy pdf 
centred on z and normalized to one on the compact subset C. The particular choice 
of a Cauchy function is suggested by a simple heuristic reasoning and confirmed by 
empirical tests. Let / be the payoff of a vanilla (call or put) option with D = 1. In 
order to compute Equation ((TJ, we should integrate a function that is of the form 
of the product of a max(-, ■) with a Gaussian pdf. As a consequence, there will be 
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a region Zn+i G {—oo, Ziow] U [z^p, +00) in which the integrand is (almost) zero and 
a region in which we expect it to be shghtly wider than a Gaussian one. We will 
come back to this point at the beginning of Section I4.2.1I 

As global benchmarks, we will quote results obtained by means of a standard Monte 
Carlo random walk (MCRW) technique. In other words, we sample N i.i.d. paths 
{Zq, Z^'-\Ti), . . . , Z^'-\Tn+i)}i=i,...,N , built up by iterating, for alH = 1 to n + 1 

Z«(Ti) = Z(')(Ti_i) + AAt + VAtS^f z = 1, . . . , AT, (17) 

where the ^ are i.i.d. standardized D-dimensional Gaussian random variables. 

Furthermore, in order to compare the PITP algorithm with a method similar in 
spirit, we implemented, in the case D = 1, a stratification-like algorithm based on Levy 
recursive construction of Brownian motion, the so called Brownian bridge. Details about 
this testing algorithm, that we will refer to as the Brownian bridge with stratification 
(BBST), are reported, for completeness, in [Appendix A[ In some cases, the algorithms 
are implemented by doubling the number of paths using antithetic variance reduction [2]. 
Whenever this is done, the corresponding algorithm is pre-fixed by the letters AT. 

4.2. Unidimensional asset 
4.2.1. Asian option 

The fair price for a discretely sampled Asian call option on an unidimensional asset 

is 

OAiZo) = E[/^(Z(T„+i),...,Zo)] 
/a(Z(T„+i),---,^) = e-'-^maxKEri)' exp{Z(T,)})/(n + 2)-ir,0} 
where K is the strike price and T the maturity. The parameters used in the numerical 
simulation are: Zq = log 100, r = 0.095, a = 0.2, t = 0, T = 1 year and n + 1 = 100. 
For simplicity, we omit the labels in the definition of Zq and a for all the unidimensional 
assets. Moreover, we consider K = 60, 100, 150 in order to test the performances of our 
algorithm when the option is ITM, ATM and OTM, respectively. 

We justify the choices made about the integration domain discussed in Section l^?11 as 
follows. Let us approximatively trace the shape of the integrand function g{zn+i)S{zn+i) 
in Equation (jllll for Asian call options. In Figure H we show the results obtained for an 
ATM, an ITM and an OTM option. Error bars correspond to one standard deviation 
Monte Carlo errors. It is particularly interesting to study the support of the integrand 
function. Actually, for ITM and ATM options, the values of Zn+i for which the function 
is considerably different from zero are more or less centred at Zo + (r — cr^/2)T, as can be 
seen from Figure H On the other hand, for OTM options, the lower bound is ^ logK. 
Hence, we can exploit this property as a rule of thumb to reduce the external integration 
to a (small) interval significantly contributing to the integral and to eventually perform 
importance sampling with an appropriate pdf. That is why we adopt particular values 
for the centre z of the integration domain. 
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Figure 1. Shape of the integrand function g{zn+i)£{zn+i) of Equation for an 
Asian cah option, showing how the support and the value of the maximum change 
when considering in-the-money (top left), at-the-money (top right) and out-of-the- 
money (bottom left) options. 
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Let us now come to the discussion of the numerical results. In Table ^ we report 
option prices and relative numerical errors as obtained by means of our path integral- 
based algorithms, as well as by the benchmarks. In the PITP and BBST cases, the 
number of integration points Uint is set to 200 and for each point we generate = 1000 
random paths. In the cases of MCRW and of pure Monte Carlo path integral with flat 
(PIFL) or Cauchy (PICH) sampling, the total number of paths is 2 x 10^, so that we 
compare results obtained with the same number of random calls and the comparison 
is meaningful. In the lower part of the table, we present the results of some of the 
algorithms improved by the implementation of the antithetic variables technique. We 
notice that all path integral prices are in very good agreement with the ones obtained 
via random walks and BBST. As a further cross-check, we compared our path integral 
predictions with the results of the method developed in |H|, finding perfect agreement. 
^From the point of view of variance reduction, the PITP algorithm turns out to be 
the method that performs best, especially when pricing ATM/OTM options. This 
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Table 1. Numerical values for an Asian call option price obtained via different 
algorithms for the parameters Sq — 100, r = 0.095, cr = 0.2, T — 1 year and n+1 = 100. 
Errors correspond to one standard deviation. 

ITA/P ATM*^ OTA/F 



Price Error Price Error Price Error 



MCRW 


40.830 


0.025 


6.899 


0.019 


0.0054 


0.0005 


BEST 


40.824 


0.018 


6.886 


0.015 


0.0058 


0.0001 


PITP 


40.811 


0.019 


6.876 


0.015 


0.0057 


0.0001 


PICK 


40.767 


0.040 


6.873 


0.019 


0.0059 


0.0001 


PIFL 


40.758 


0.105 


6.880 


0.026 


0.0057 


0.0001 


AT-MCRW 


40.836 


0.002 


6.909 


0.008 


0.0053 


0.0003 


AT-PITP 


40.832 


0.004 


6.901 


0.004 


0.0060 


0.0001 


AT-PICH 


40.775 


0.031 


6.878 


0.008 


0.0058 


0.0001 


^ In-the-money, K — 


60. 










^ At-the-money, K ^ 


: 100. 










Out-of-the-money, K = 150. 











means that, when the integrand is non-zero only in a region far from Zq + {r — (T^/2)T, 
the standard MCRW generates many paths that are not relevant for the mathematical 
expectation. Furthermore, the PITP and the BBST techniques give essentially the same 
results, thus confirming our suspicion that the strategy of fixing the end point before 
generating paths plays the crucial role in the variance reduction. Let us stress that the 
flat integration gives a worse variance and that PICH and PIFL algorithms perform 
best only out of the money. 

To conclude, let us comment about the estimate of the numerical error for the PITP 
and the BBST algorithms. Errors in Equation (jl5|) result from the combination of the 
Monte Carlo errors on each end point, -Zn-i-i- To estimate the error associated with the 
finiteness of Uint required by the deterministic integration, we analyzed the stability of 
the price with respect to the number of integration points. In Figure El we show the 
prices obtained according to this procedure. It can be seen that the fluctuations of the 
price value due to the choice of Uint become negligible with respect to the width of the 
error bars, as the value of Uint increases. This is why we consider the relevant source 
of numerical error as related to the Monte Carlo part of the integration and we set 
Hint = 200 in our simulations. 

4-2.2. Up- out barrier options 

In this section we consider barrier options of European style, i.e. whose exercise is 
possible only at the maturity. In particular, we price so called knock-out up options. 
The payoff at maturity Tn+i = T is a functional of the whole path {5(5); < s < T} 
and has the following expression 

f^[S] = e""^ max(^(T) - K, 0)I,>t + e-"("-*)i?I,<r, (19) 
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Figure 2. Prices and error bars for an at-thc- money Asian call option as a function 
of the number of points for external integration riint. As ni„t increases, the error bars 
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where U is the upper barrier, r = inf{s > : S{s) > U), R is a fixed cash rebate and 
is the characteristic function of the set A. In our simulation we set R = 0. It is 
known j21I that, whenever we discretize this continuous time problem, setting 

n+l 

fuiZr, ...,Zo)= e-^^ max(e^(^) - ^, 0) J] I^mXiogC/, (20) 

1=0 

we overestimate the price of up-out options. Actually, we do not take into account the 
possibility that the asset price could have crossed the barrier for some t G]Tj, Tj+i[, i = 
to n. A strategy to obtain a better approximation is the following. We check, at each 
time step Tj = iAt, and for all Monte Carlo paths, whether exp{Z^^)(Tj)} has reached 
the barrier U. If not, we compute the value 

-^(logf/-ZW(T._0)(logf/-^«(T.))j , (21) 

and we extract a random variable from a Bernoulli distribution with parameter pi. If 
the result is 1, the barrier value has been reached in the interval ]Tj_i, Tj[ and the price 
associated to that given path is zero. Otherwise, the simulation is carried on. This 
technique is widely discussed in [211 123 1211 ■ 

In Table El we report prices, with the corresponding numerical errors. We price an 
ATM option, K = 100, and an OTM option, K = 130, with U = 150 and U = 200, all 
the other parameters as in Section 14.2.11 It is particularly evident that the precision of 
PITP exceeds that of MCRW, both for ATM and OTM options. 

4-2.3. Reverse cliquet options 

The last one- dimensional example we consider is represented by the reverse cliquet 
option, whose payoff function is 



fRc{Z{Tn+i), . . . , zo) = e '■'^ max 



.(22) 

S{Ti)=exp{Z(TO} 
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Table 2. Numerical values for the price of up-out barrier call options obtained through 
different algorithms for the parameters Sq = 100, r = 0.095, a ~ 0.2, T = 1 year and 
n + 1 = 100. Errors correspond to one standard deviation. 

K = 100 K = 130 



U = 150 U = 200 U = 150 U = 200 

Price Error Price Error Price Error Price Error 

AT-MCRW 9.087 0.012 12.853 0.015 0.647 0.004 2.353 0.011 

AT-PITP 9.088 0.008 12.830 0.001 0.638 0.002 2.336 0.001 

AT-PICH 9.099 0.016 12.815 0.014 0.647 0.002 2.333 0.003 



Table 3. Reverse cliquet option fair price for 5*0 = 100, r = 0.09, a ~ 0.3 and 
T/n = 1/12 year. Errors are not quoted in 





n = 


4 


n = 


12 


n = 


24 


n = 


36 




Price 


Error 


Price 


Error 


Price 


Error 


Price 


Error 


AT-MCRW 

AT-PITP 

AT-PICH 

m 


0.0574 
0.0572 
0.0572 
0.0574 


0.0001 
0.0001 
0.0001 


0.1223 
0.1225 
0.1218 
0.1222 


0.0001 
0.0002 
0.0002 


0.1993 
0.1992 
0.1990 
0.1990 


0.0002 
0.0003 
0.0002 


0.2611 
0.2611 
0.2606 
0.2609 


0.0002 
0.0003 
0.0003 



The option is characterized by the number of periods n, the minimum amount (the 
floor) F, and the maximum payable coupon (the cap) C. This option is especially 
valuable when positive performances are more probable. 

We have tested our algorithms by choosing the parameters' values as in 0, whose 
numerical results for option prices are reported in Table El for the sake of comparison. 
The floor F has been fixed equal to zero, while the cap C is set equal to no, with 
c = 0.04. The spot is 5*0 = 100, r = 0.09 and a = 0.3. We also change the maturity T 
by fixing T/n = 1/12 year and letting n = 4, 12, 24, 36. In Table 01 we show the values 
corresponding to different values of n. Integration on 2;„_|_i for the AT-PITP algorithm 
is performed on a symmetric interval of width 8cr\/T centred on Zq + (r — cr^/2)T. 
Once again we observe accurate path integral pricing and a good agreement between 
the results of the various algorithms. 

4.3. Multidimensional assets 

In this section we study the performances of path integral pricing in the case of options 
on multidimensional assets S = {Si, S2, ■ ■ ■ , Sd). As an example, we price an Asian call 
option on the basket X whose value at time t is obtained by linearly combining the 
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values of the components of S: 

EZlC.^ = 1 (23) 
ai>0 \fi = 1 to D. 

Consequently, we have 

n+l D \ 

—— 5^ 5^ a. exp{Z\T,)} - K,o\ .(24) 

i=0 1=1 J 

In order to compare the multidimensional performance of all the algorithms 
introduced in Section El we need to choose D such that it still makes sense to perform 
a deterministic integration over Z(T„+i) = log5'(T). However, we expect a gain in 
competitiveness of pure Monte Carlo integration (PIFL, PICH), as the deterministic 
approach gradually loses its attractive features as the dimension increases. That is why 
we choose a three-dimensional correlated asset. 

All the tests are performed by setting r = 0.095 and considering a maturity of 
T = 1 year with a time discretization of 100 time steps (n + 1 = 100). Unless otherwise 
specified, pij = 0.6 for any i ^ j, and cr^ = 0.2 for all = 1 to D = 3. Path 
integral integration is limited to the compact subset C described in Section 14.11 In the 
special case of PITP, we consider 1000 Monte Carlo samples for each end point. It is 
worth noticing that, if we choose ni integration values for each dimension, the total 
number of required deterministic integration points grows exponentially as nf. Thus, 
an apparently poor-quality unidimensional integration with ui = 10 consists indeed in 
evaluating the integrand function on 10^ points. 

The first test concerns the convergence of the deterministic integration. We set 
K = 120, So = (100, 90, 105) and we compare prices obtained with ui = 6, 8, 10, i.e. with 
216 ■ 10^, 512 • 10^, 10^ total integration points for Zn+i- In Table ID we show the prices 
thus obtained with their one standard deviation errors, together with the ratio between 
one standard deviation Monte Carlo error and the price corresponding to ni = 10 
(fourth column), and the percentage difference between Price(ni) and Price(?T,i = 10) 
(fifth column). As in the unidimensional case, changes in prices due to the number of 
deterministic integration points fall well inside the Monte Carlo 95% confidence interval 
[Price — 1.96 Error, Price + 1.96 Error]. 

As a second test, we studied the drawbacks of performing pure Monte Carlo (PIFL, 
PICH, MCRW) simulations with a small number of Monte Carlo paths (10^) finding 
that in these cases MCRW works best, since the path integral fails to efficiently explore 
the support of the integrand function. Results corresponding to pij = 0.8 for any i j, 
K = 110, (Tfc = 0.2 + 0.02 ■ (A; — 1) and obtained with two different choices for the centre 
z of the integration hyper-cube C and two different spot 5*0 are reported in Table El The 
case 5*0 = (100, 95, 80) corresponds to an OTM option, while when 5*0 = (107, 109, 114), 
corresponds to an ATM option. Convergence results are poor as the estimated price 
depends on the integration region and the Monte Carlo error is large. 
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Table 4. Prices, one standard deviation Monte Carlo errors, percentage errors 
and prices' percentage differences of an Asian basket call option, according to the 
multidimensional PITP algorithm, with K = 120. The reference for prices' percentage 
differences is Price(ni = 10). 



ni 


Price 


Error 


Relative 


Percentage Price 








Error 


Difference 


10 


0.306 


0.003 


0.99% 




8 


0.310 


0.005 


1.48% 


1.23% 


6 


0.323 


0.008 


2.38% 


5.9% 



Table 5. Prices and errors for Asian basket call options according to the algorithms 
PIFL, PICH, MCRW with a small set (10'') of Monte Carlo samples. Errors correspond 
one standard deviation. K — 110 





ATM {So = (107, 109, 114)) 


OTM {So = 


(100,95,80)) 




e^=(105,110,115) 


e*=(100, 120,115) 


c*=(110,110,110) 


e^=(120,120,110) 




Price Error 


Price Error 


Price Error 


Price Error 


MCRW 

PIFL 

PICH 


7.5 0.1 
6.3 0.5 
6.8 0.4 


7.5 0.1 
7.8 0.6 
8.1 0.4 


0.83 0.03 
0.77 0.09 
0.80 0.09 


0.83 0.03 
0.78 0.08 
0.56 0.07 



Once we take care of choosing a sufficient number of Monte Carlo paths ^, we 
compare the performance of our algorithms in the cases of ATM and OTM options. 
We report the results of such analysis in Table where 5*0 = (100, 90, 105) and the 
strike K is varied to perform ITM/ATM {K = 100) and OTM {K = 140) pricing. We 
report results obtained with two differently centred integration intervals, in order to 
show that the chosen number of samples is enough to guarantee stability of integration 
to (relatively small) changes of the integration hypercubeff. When compared to Tabled 
the results shown in Table El present some similarities and some differences. As in the 
unidimensional case, the path integral is still a valuable choice to price OTM options, 
prices being in agreement with the benchmark MCRW and errors smaller, especially 
with a Cauchy pdf sampling (PICH). As expected, however, the external deterministic 
integration (PITP) has effectively lost its attractive properties, PIFL and PICH giving 
more precise confidence intervals. Let us stress that, when the dimensionality increases, 
the path integral performance for ATM options worsen. Even if we perform tests 
with 5 ■ 10^ and 10® Monte Carlo samples, whenever we price ITM/ATM options, the 
path integral method is a bit less precise than the random walk, the central value 
being nevertheless compatible with the benchmark results. Because the performances 

Just take as reference the deterministic integration with the rule of thumb of setting at least 6 
integration points for each dimension and 1000 MC paths for each end point, i.e. 6^ ■ 10'^ total samples, 
ff It is important to recall that, if we perform integration on an interval whose spot values are too low, 
we will have an under-estimation of the price. 
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Table 6. Prices and errors for Asian basket call options according to the algorithms 
PITP, PIFL, PICH, MCRW with ni = 6 and 216000 total Monte Carlo paths. 
So = (100,90, 105). Errors correspond to one standard deviation. 







ATM {K 


= 100) 






OTM {K 


= 140) 




e^=(110,100,110) 


e*=(100, 100,100) 


e*=(140,140,140) 


e*=(130, 130,130) 


Price 


Error 


Price 


Error 


Price 


Error 


Price 


Error 


MCRW 


5.29 


0.02 


5.29 


0.02 


0.0049 


0.0004 


0.0049 


0.0004 


PITP 


5.33 


0.04 


5.28 


0.04 


0.0051 


0.0003 


0.0049 


0.0003 


PIFL 


5.37 


0.06 


5.41 


0.07 


0.0048 


0.0003 


0.0048 


0.0002 


PICH 


5.26 


0.03 


5.28 


0.03 


0.0048 


0.0001 


0.0050 


0.0001 



for OTM options are good, we infer that, to increase the accuracy in pricing high 
dimensional ITM/ATM options, we should consider an hyper-cube C wider than four 
standard deviation. 

5. Greeks 

In the present section, we report our results concerning the computation of the Greeks 
for unidimensional Asian and barrier knock out options. Actually, it is interesting to 
compare the numerically estimated exotic Greeks with those implied by the Black and 
Scholes model for plain vanilla call options. The values of the parameters are K = 100, 
r = 0.095, a = 0.2, T = 1 year and n + 1 = 100, while for the barrier we choose U = 150. 

In Figure El we show the price of the option and the Greeks delta, gamma, 
vega and theta as functions of the spot price S. In our approach the Greeks are 
numerically computed using a finite difference method applied to the option values 
returned by the path integral with trapezoidal integration. The error bars are obtained 
by propagating the (one standard deviation) numerical error of the path integral option 
prices. Therefore, the values of the Greeks are more accurate than those obtained with 
other approaches when the path integral pricing outperforms existing techniques. 

As expected, the qualitative behaviour of prices and Greeks for Asian call options 
does not differ significantly from that of plain vanilla ones. The shift in prices is due 
to the fact that in the Asian payoff the role of S{T) is played by the mean value of 
S along the path. A lower price is therefore a straightforward consequence of the fact 
that we have E[S'(Ti) + • • ■ S{Tn+i)]/{n + 1) < E[S'(T„+i)]. The Greeks do not coincide 
exactly, but the relevant features, such as the sign of the derivative, are preserved. The 
situation is completely different for the barrier options. First of all, we expect that for 
small values of S and with our choice for the max barrier value, U = 150, the results 
overlap the European ones. This can be easily verified from Figure El and considered 
as a check of the consistency of our numerical results. The profile of the price graph 
is characterized by changes both of the monotonic properties and of the concavity, as 
shown in FigureOl This reflects in the change of sign of delta and gamma. The behaviour 
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Figure 3. Option price and Greeks for plain vanilla (- 
knock out (O ) call options. 
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of vega can be explained by noticing that, for S U, an increase of the a means an 
increase in the width of the distribution of S{T), that reaches higher values without 
reaching the barrier, so 80/ da > 0. Conversely, when S and a grow, S{T) reaches 
the barrier more frequently and the option loses value. Analogous reasoning applies to 
theta, i.e. —dO/dT, where the role of a is played by the maturity T. However, in this 
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case, the presence of the minus sign in the definition imphes < for S <^ U and 
O > otherwise. 

6. Conclusions and outlook 

In this paper we have shown how the path integral approach to stochastic processes can 
be successfully applied to the problem of pricing exotic derivative contracts. Numerical 
results for the fair price and the Greeks of a variety of options have been presented 
and compared with those obtained by means of other standard (such as the Monte 
Carlo Random Walk) and non standard (see 0) approaches employed in quantitative 
finance. With respect to the original formulation of [0] , the method has been generalized 
in order to cope with options depending on multiple and correlated underlying assets. 
Concerning options depending on a single asset, it has been shown that the algorithm can 
provide very precise results, especially when pricing ATM and OTM options. This is due 
to an appropriate separation of the integrals entering the path integral pricing formula 
and, more importantly, to a careful simulation of random paths with fixed end points, 
able to probe the regions contributing to the dominant part of the payoff functions. As 
far as (multidimensional) basket options are concerned, while the standard Monte Carlo 
simulation turns out to be more efficient for ITM/ATM options, our approach exhibits 
better performances for OTM options. 

In all the cases, the computational time is essentially the same required by a 
standard Monte Carlo calculation. The algorithm is general and could be extended 
to price other types of exotic contracts. 

As a future important development, it would be interesting to explore the feasibility 
of an application of our framework to derivative pricing approaches beyond Black-Scholes 
and dealing with non-Gaussian features of financial time series, such as Bouchaud- 
Sornette residual risk minimization 12H1 122] , stochastic volatility models ISOIISIIE^I, 
multi fractal random walks [SSIISII, GARCH processes |3SI and other generalizations of 
the log- normal dynamics [^HIEZIEHI- 
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Appendix A. Stratification and Brownian bridge 

We describe here the algorithm used to test our hints about the reasons of the 
good performances of path integral with deterministic integration when pricing path 
dependent options on unidimensional assets. In order to improve the numerical 
precision, it is necessary to lead Monte Carlo paths to a region in which the payoff 
function is different from zero. This is the advantage of performing the external integral 
in Equation (fTTjl in a clever way and to drive paths toward some fixed end points. The 
algorithm here described exploits this idea by means of a backward construction of the 
underlying process (Brownian bridge), instead of using a path integral method. For 
simplicity, since this method has only been used in the case D = 1, we outline here the 
construction of a unidimensional Brownian bridge only. 

The aim is to describe the law of Z{t), with t e [0, t'], once Zq and Z{t') are known. 
By setting 

Z{t) = Zo+ ^^^'\;^' t + e{t), (A.l) 

it is easy to see that e{t) is Gaussian with zero mean and variance t{t' — t)a'^/t'. It 
is worth noticing that, while LaM[Z{t)\Z(t')] 7^LawZ(t), e{t) is independent of Z{t'). 
Hence, given Zq and Z(T), we construct Monte Carlo paths Zq, Z(Ti), . . . , Z(T) by 
recursively applying Equation ()A.1|) from t' = T, t = down to t' = T2, t = Ti. As 
for the path integral framework of Equation p3|) and the standard Euler discretization 
scheme of Equation (|T7jl . sampling a path is equivalent to generating i.i.d. standardized 
Gaussian variables. The main difference with the path integral method is that here, once 
we have the l-th. Gaussian sample {^'■'^(Ti), . . . ,£:(')(T„)}, we are nevertheless forced to 
simulate the path recursively, which means that we cannot infer Z^'-^Ti) if we have not 
sampled Z*^')(Tj+i). In the path integral case, instead, extraction is straightforward by 
means of Equation (fT^ . 

Let us now explain how to compute Equation by means of the Brownian bridge. 
By setting /(zn+i, ?/„..., yi, zq) = f{zn+i, z„ . . . , zq) with the replacements 

Zi = Zq^ Ti + yi, 1 = 1, ...,n, 

we have 

E[/(Z(T„,+i), ...,Zq)]= E[/(Z(T„+i), £(T„), . . . , £(Ti), ^0)] 



d2n+lPn+l(^n+l) 



« n 

/ T\idyiPe,iiyi))fi^n+l,yn,---,yi,ZQ) (A.2) 



dZn+lpn+liZn+l)'E.[fiZn+l,e{Tn), e(Ti), Zq)], 

where pn+i is the pdf of Z(T) and p^^j is the pdf of e(Ti), for i = 1 to n. 

It is now straightforward to see that Equation ()A.2j) has the same form as 
Equation (|TT|). Therefore, we can approximate integration over Zn+i as in Equation (fT^ 
and evaluate via Monte Carlo the inner mathematical expectation. 
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By means of Equation (jl4p . we are, in some sense, stratifying the domain V of 
the random variable Z(T„+i) by dividing it into disjoint sub-sets T>i (here the sub-sets 
reduce to the points G M of Equation (fT^ ). For each sub-set, we then compute the 
inner integral by forcing Z{Tn+i) & T^i. It is possible to show |23| that this procedure 
may lead to variance reduction. This way of proceeding has the same qualities and the 
same limitations as the PITP, yielding less accurate results for multidimensional assets. 
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